First-principles study of the polar O— terminated ZnO surface 
in thermodynamic equilibrium with oxygen and hydrogen 
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Using density-functional theory in combination with a thermodynamic formalism we calculate 
the relative stability of various structural models of the polar O-terminated (0001)-O surface of 
ZnO. Model surfaces with different concentrations of oxygen vacancies and hydrogen adatoms are 
considered. Assuming that the surfaces are in thermodynamic equilibrium with an O2 and H2 gas 
phase we determine a phase diagram of the lowest-energy surface structures. For a wide range of 
temperatures and pressures we find that hydrogen will be adsorbed at the surface, preferentially 
with a coverage of 1/2 monolayer. At high temperatures and low pressures the hydrogen can be 
removed and a structure with 1/4 of the surface oxygen atoms missing becomes the most stable 
one. The clean, defect-free surface can only exist in an oxygen-rich environment with a very low 
hydrogen partial pressure. However, since we find that the dissociative adsorption of molecular 
hydrogen and water (if also the Zn-terminated surface is present) is energetically very preferable, it 
is very unlikely that a clean, defect-free (0001)-O surface can be observed in experiment. 

PACS numbers: 68.47.Gh, 68.35.Md, 68.35.Bs, 71. 15. Mb, 82.65.+r 
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I. INTRODUCTION 

To understand the structure, composition and stabil- 
ity of polar surfaces on a solid theoretical basis is one of 
the challenging problems in surface sciencei The most 
interesting polar surfaces are so called "Tasker-type(3)" 
surfaces 2 which are formed by alternating layers of op- 
positely charged ions. Assuming a purely ionic model 3 
in which all ions are in their formal bulk oxidation state, 
such a stacking sequence creates a dipole moment perpen- 
dicular to the surfaces which diverges with slab thickness, 
and with simple electrostatic arguments it can be shown 
that the surface energy will diverge with sample size^ To 
quench the dipole moment and to make the polar surfaces 
stable, a redistribution of charges in the surface layers has 
to take place.- For various polar surfaces different mecha- 
nisms to accomplish the charge compensation have been 
observed^ however, in many cases the underlying sta- 
bilization mechanism is very controversially discussed in 
the literature. 

One of the most widely investigated examples of 
Tasker-type(3) polar surfaces are the two basal planes of 
ZnO: the O-terminated (OOOT)-O and the Zn-terminated 
(OOOl)-Zn surface. The two surfaces are the terminating 
planes of a stacking sequence of hexagonal Zn and O 
layers along the crystallographic c-axis with alternating 
distances of i?i=0.6lA and i?2=l-99A. In this case, as 
can be easily shownii, the polar ZnO surfaces are only 
stable if the O-terminated face is less negative and the 
Zn-terminated surface layer less positively charged com- 
pared to the formal bulk oxidation state by a factor of 
Ri/(R\ + i?2)~l/4. In principle, three different scenar- 
ios are conceivable to accomplish this charge redistribu- 
tion: The ionic charge of the surface ions may be reduced 
from ±2 to ±3/2, which may be regarded as an "electron 
transfer" from the O- to the Zn-terminated surface (la). 
As a result, partially occupied surface bands will appear 



with a 3/4 filled 0-2p band at the (0001)-O and a 1/4 
filled Zn-4s band at the (OOOl)-Zn surface. This is of- 
ten referred to as "intrinsic surface state compensation" 5 
or as "metalization of the polar surfaces" £ However, 
whether a true metallic state is present will depend on 
the dispersion of the partially occupied bands. Addition- 
ally, in a second step, the surface may reconstruct and 
undergo a distortion in which, for example for the O- 
tcrminated surface, four surface atoms combine in such 
a way that an unoccupied 2p-band splits from the other 
eleven occupied 2p-bands and the surface becomes in- 
sulating again (lb). Secondly, the charge reduction of 
the surface layers may take place by removing 1/4 of 
the surface ions (II). The so created vacancies may be 
ordered and form a reconstruction or may be randomly 
distributed. Finally, charged species may be adsorbed to 
reduce the formal oxidation state of the surface ions (III) . 
For example, water may dissociate and protons (H + ) and 
hydroxyl groups (OH - ) adsorb on every second O and 
Zn surface ion, respectively^ All three scenarios repre- 
sent idealizations of the charge compensation process. In 
general, any combinations of the three mechanisms are 
conceivable, like the simultaneous formation of vacancies 
and partially filled bands, as long as the charge com- 
pensation rule is obeyed. The surface structure finally 
realized for a specific polar surface will then be the one 
with the lowest surface energy. 

For ZnO it was believed for a long time that both 
polar surfaces exist in an unreconstructed, truncated- 
bulk-like state. After standard preparation proce- 
dures both surfaces show regular (lxl) pattern in low- 
energy electron diffraction (LEED) 10 and other diffrac- 



tion experiments 
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Some evidence for missing Zn 



ions on the (OOOl)-Zn surface was found in gracing in- 
cidence X-ray diffraction (GIXD) 11 , however, for the 
(0001)-O surface no evidence for substantial amounts of 
surface oxygen vacancies was detected in GIXD 1 ! and 
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low-energy alkali- ion scattering (LEIS)^ 

For ideal, truncated-bulk-likc surface terminations 
only mechanism (la) can explain the stability of the po- 
lar ZnO surfaces. Consequently, in most theoretical first- 
principles studies of the polar ZnO surfaces ideal surface 
terminations together with partially filled surface bands 
were assumedi&ii Studies exploring the other two stabi- 
lizations mechanisms are very scarce. In a pioneering ab- 
initio study Wander and Harrisoni4 investigated whether 
the polar surfaces may be stabilized by the dissociation of 
water and the adsorption of H + and OH~ groups accord- 
ing to mechanism (III). They found this energetically un- 
favorable compared to situation (la). However, instead of 
1/2 monolayers, which would be the ideal configurations 
for band filling, only full monolayers of H + and OH~ were 
considered, thereby overcompensating the needed charge 
transfer between the polar surfaces. In addition, only one 
specific adsorption site for the H + and OH~ groups was 
probed. 

Meanwhile two recent experiments have created con- 
siderable doubt that the polar ZnO surfaces really exist 
in a clean, unreconstructed state. With scanning tun- 
neling microscopy (STM) it was showni^ii& that the Zn- 
terminated surface is characterized by the presence of 
nanoscaled, triangular islands with a height of one ZnO 
double-layer. The shape of the islands is size-dependent 
and typical pit structures appear for larger islands. Since 
the step edges are O-terminated, the high step concen- 
tration leads to a significant decrease of Zn ions in the 
surface. A rough analysis of the island and pit size dis- 
tribution yielded that approximately 1/4 of the Zn ions 
is missing, in agreement with mechanism (II). With de- 
tailed density-functional theory (DFT) calculationsi^ii 
it was confirmed that a crystal termination with trian- 
gular shaped islands and pits is indeed lower in energy 
than the perfect bulk-truncated surface for a wide range 
of oxygen and hydrogen chemical potentials. Under H 
rich conditions, structures with up to 1/2 monolayer of 
hydroxyl groups were even more stable, indicating that 
the actual surface morphology will sensitively depend on 
the chemical environment. 

On the other hand, for the O-terminated polar sur- 
face no such island and pit structure was found with 
STMii& However, with He scattering (HAS) it was 
discovered^ that at ultrahigh vacuum (UHV) condi- 
tions O-terminated surfaces with (lxl) LEED and HAS 
diffraction pattern are always hydrogen covered, whereas 
after a careful removal of the hydrogen a (1x3) structure 
is found. The (1x3) spots are best visible in HAS, but 
under certain conditions can also be observed in LEED. 18 
The H-free surface is very reactive and dissociates molec- 
ular hydrogen and water and therefore exits only for a 
limited time even at UHV conditions. In a subsequent 
stud y 19 i 20 i 21 i 22 CO was used as a probe molecule to dis- 
tinguish between clean and hydrogen saturated surfaces. 
By comparing calculated CO adsorption energies for dif- 
ferent surface structures with experimental results it was 
confirmed that the polar O-terminated surface is usu- 



ally hydrogen covered whereas a clean (lxl) (0001)-O 
surface is very unlikely to exist. 

In previous studies the H coverage of the polar O- 
tcrminated surface was not observed most likely because 
LEED and X-rays are not sensitive to hydrogen. How- 
ever, it is not clear how the structures found in the HAS 
study, Ref. HI lead to a stabilization of the polar O- 
terminated surface. A full hydrogen monolayer would 
overcompensate the charge transfer so that again par- 
tially occupied bands have to be present, and the nature 
of the H-free (1x3) structure is still unknown. There 
is some evidence from X-ray photoelectron spectroscopy 
(XPS) 18 that oxygen vacancies are involved, but 1/3 of 
the oxygens missing is also not the expected vacancy con- 
centration to stabilize the surface. 

Motivated by these experimental findings we explore 
in the present paper the competition between the three 
stabilization mechanisms in a very general way. The 
main focus will be on the O-terminated (0001)-O sur- 
face, and we take an approach very similar in spirit to 
the investigation of the Zn-terminated surface in Ref. ITEI 
and For a series of surface models we determine the 
total energies and the fully relaxed atomic structures us- 
ing a first-principles DFT approach. Surface structures 
with various oxygen vacancy concentrations and different 
amounts of adsorbed hydrogen are considered, including 
structures corresponding to the three ideal stabilization 
scenarios and structures compatible with the HAS obser- 
vations. 

Static total-energy DFT calculations only give results 
for zero temperature, zero pressure and for surfaces in 
contact with vacuum. However, the actual lowest-energy 
structure of the (0001)-O surface will depend on the en- 
vironment and can change with temperature T, pressure 
p and exposure to O2 and H 2 gas phases. Therefore, to 
determine the equilibrium structure and composition of 
the surface at finite temperature and oxygen and hydro- 
gen partial pressures, we combine our DFT results with 
a thermodynamic description of the surfaces. To take de- 
viations in surface composition and the presence of gas 
phases into account, we introduce appropriate chemical 
potentials^ and calculate an approximation of the Gibbs 
free surface energy^ Depending on the chemical poten- 
tials we then determine the surface structure with the 
lowest free energy which allows us to construct a phase 
diagram for the surface. If we assume that the surface 
is in thermodynamic equilibrium with the gas phases, we 
can relate the chemical potentials to a given temperature 
T and pressure p. In this way we are able to extend our 
zero temperature and zero pressure DFT results to ex- 
perimentally relevant environments, thereby bridging the 
gap between UHV-like conditions and temperatures and 
gas phase pressures that are typically applied, for exam- 
ple, in catalytic processes like the methanol synthesis^ 
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II. COMPUTATIONAL APPROACH 
A. Thermodynamics 

In this section we will give a brief description of 
the thermodynamic formalism which we have used 
to determine the most stable structures of the po- 
lar O-terminated ZnO surface. The formalism has 
been successfully applied in several previous surface 
studie3i£*2iii2§i2L22iSi and is described in more detail in 
Refs.0and|2i 

The general expression for the free energy of a surface 
in equilibrium with particle reservoirs at the temperature 
T and pressure p is given by^S 

j(T,p) = j(G(T,p,{N i })-Y^N i ^ i (T,p) ] j , (1) 

where G(T,p, {Ni}) is the Gibbs free energy of the solid 
with the surface of interest, A is the surface area, and 
Ni are the chemical potentials and particle numbers of 
the various species. In contrast to the usual convention in 
macroscopic thermodynamics we define here the chemical 
potentials per atom rather than per mole. For the study 
of the polar O-terminated ZnO surface in contact with 
an oxygen and a hydrogen gas phase we have to consider 
the three chemical species i = Zn, O and H. 

For simplicity we have assumed two independent reser- 
voirs for O2 and H2 with a common pressure p. Experi- 
mentally, it is more likely that a mixture of O2 and H2 is 
present. In this case, the pressure p in Eq. has to be 
replaced by appropriate partial pressures po 2 and pn 2 - 
However, in the present study we will keep the restric- 
tion of separate reservoirs in the sense that we do not 
allow O2 and H2 to react to H2O, which would be the 
case in full thermodynamic equilibrium. This is justified 
by arguing that the reaction barrier for the formation of 
H2O is high enough that the reaction plays no role on 
the time scales of interest. 

In thermodynamic equilibrium the chemical potentials 
would be determined uniquely by the temperature T, the 
pressure p and the total particle numbers of the solid and 
the gas phases. The surface structure, here represented 
by the particle numbers Ni, would then be determined 
by a unconstrained minimization of the surface free en- 
ergy, Eq. (I). However, this is not very practical to do. 
Therefore we will take a different approach: We calcu- 
late the surface free energy of a series of model surfaces 
with different structures and compositions as a function 
of the chemical potentials. For given chemical potentials 
we predict which surface structure is the most stable one 
by searching for the surface model with the lowest surface 
free energy. In a second step, the chemical potentials are 
then related to actual temperature and pressure condi- 
tions by assuming that the surface is in thermodynamic 
equilibrium with the gas phases. 

In our calculations all surfaces are represented by pe- 
riodically repeated slabs so that the Gibbs free energy 



G(T,p, {Ni}) refers to the content of one supercell. Since 
ZnO is not centrosymmetric, slabs representing the po- 
lar ZnO surfaces are inevitably O-terminated on one side 
and Zn-terminated on the other side. It is therefore not 
possible to assign unique surface energies to the two po- 
lar surface terminations. Only the sum of the surface 
energies and thereby the cleavage energy are well defined 
quantities. However, in the present study we are only in- 
terested in the relative stability of O-terminated surfaces 
with different structures and compositions. The Zn-face 
of the slabs is unchanged in all calculations. Therefore 
we relate all energies relative to a reference state which 
we have taken to be the ideal, truncated-bulk termina- 
tion. We define the change of the cleavage energy A7 as 
the difference of Eq. JIJ for a slab with a chosen surface 
structure and a slab with ideal surface terminations: 

A7(7» - j(G s Zl(T,p,N v ,Nn)-G^ b (T,p) 

+ N vf i {T,p)-Nnm(T,p)) . (2) 

Here, and G^ b are the Gibbs free energies of the 

supercells with the model surface and the reference con- 
figuration, respectively, and A is now the area of the 
surface unit cell. Since we only consider structures of the 
O-terminated surface with O-vacancies and adsorbed H 
atoms, only the number of O-vacancies JVy (= difference 
of the number of O-atoms in the slab with the model sur- 
face and the reference state) and the number of adsorbed 
H atoms N-g. appears in Eq. (J3J), i.e. the chemical poten- 
tial of Zn drops out. The difference A7 is negative if a 
model surface is more stable than the ideal, truncated- 
bulk-like surface termination and positive otherwise. 

In principle we now have to calculate the Gibbs free 
energy of all slabs representing our surface models, in- 
cluding the contributions coming from changes in volume 
and in entropy. The entropy term may be calculated, 
for example, by evaluating the vibrational spectra in a 
quasiharmonic approximation^, but in practice this is 
computationally very demanding. As is apparent from 
Eq. P|l. only the difference of the Gibbs free energy of 
two surface structures enters the expression for A7. In 
Ref.|29|it was shown that the vibrational contributions to 
the entropy usually cancel to a large extent and that the 
influence of volume changes are even smaller. Therefore 
we will neglect all entropy and volume effects. The Gibbs 
free energies then reduce to the internal energies of the 
slabs and we can replace G^ and G r s fl b in Eq. © by 
the energies as directly obtained from total-energy (e.g. 
DFT) calculations. 

Finally we have to determine meaningful ranges in 
which we can vary the chemical potentials. First, there 
are upper bounds for all three chemical potentials no, 
/in and ^tzn, beyond which molecular oxygen and molec- 
ular hydrogen would condensate and metallic Zn would 
crystallize at the surface. These bounds are given by the 
total energy of the isolated molecules Eq 2 , Eg 2 and of 
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bulk Zn E^ lk (neglecting volume and entropy effects): 

Mo < 2 E ° 2 ' ^ H ~ 2 Eii2 ' ^ Zn ~ ^n' k • (3) 

In the following we will use these upper bounds as zero 
point of energy and relate the chemical potentials relative 
to the total energies of the isolated molecules: 

A/x = Mo - \ E o 2 , A^ H = MH - ^E H2 . (4) 

Furthermore we impose that the surface is always in equi- 
librium with the ZnO bulk phase. Then the sum of no 
and /izn has to be equal to the total energy E Zn o of bulk 
ZnO. Thus only one of the two chemical potentials fio 
and fiin is independent, and together with Eq. ||2J we in- 
troduce lower bounds for the chemical potentials. Using 
the energy of formation Ef of ZnO: 

El = EznO — Ezn — 2^°2 (5) 

the allowed range for the chemical potential A/io is given 
by: 

E { < Ano < . (6) 

If we assume that the surfaces are in thermodynamic 
equilibrium with the gas phases we can relate the chem- 
ical potentials A/io and A/xh to a given temperature T 
and partial pressures po 2 and pn 2 - F° r ideal gases we 
can use the well-known thermodynamic expressions^ 

A Mo (T,po 2 ) = l(vo 2 (T,P°) + k B T]n(po 2 /p )) (7) 
and 

A/in {T, ph 2 ) = l(m 2 (T,p°) + fc B Tln(pH 2 /p )) (8) 

in which p° is the pressure of a reference state and 
the temperature dependence of the chemical potentials 
ft,o 2 (T,p°) and jdff 2 (T,p ) is tabulated in thermochem- 
ical reference tables^ However, it should be noted, as 
pointed out by Finnish, that the equilibrium with the 
gas phase need not to be perfect. It is sufficient if the 
surface is in equilibrium with the bulk phase. In this 
case, the chemical potential is related to defect concen- 
trations of the bulk. For example, if oxygen vacancies 
are the dominant defects we haveM^^ 

Ano(T,p) = Mo - fc B rin(c v /co) (9) 

with the vacancy concentration cy and the oxygen occu- 
pancy of the O-lattice site of cq. 

B. Ab-initio Calculations 

The total energies of slabs representing different model 
surfaces as well as the bulk and molecular reference en- 
ergies were calculated within the framework of density- 
functional theory (DFT)M> Exchange and correlation ef- 
fects were included in the generalized-gradient approxi- 
mation (GGA) using the functional of Perdew, Becke and 



H 2 2 Bulk Zn Bulk ZnO 
Ef BE [eV] 4.50 6.38 1.12 2.84 
# t cxp [eV] 4.52 5.17 1.35 3.50 

TABLE I: Calculated binding energies Ef BE for the isolated 
H2 and O2 molecules and for the bulk phases of Zn and ZnO. 
Zero-point vibrations are not included. The experimental 
heat of formations Hi ' xp are for T=298K and p=lbar and 
are taken from Ref. 144 



Ernzerhof (PBE). 38 Normconserving pseudopotentials 3 8 
were employed together with a mixed-basis set consist- 
ing of plane waves and non-overlapping localized orbitals 
for the 0-2p and the Zn-3<i electrons^ A plane-wave 
cut-off energy of 20 Ry was sufficient to get well con- 
verged results. Monkhorst-Pack k-point meshes 41 with a 
density of at least (6x6x6) points in the primitive ZnO 
unit cell were chosen. A dipole correction 4 ^ 3 , to the 
electrostatic potential was included in the calculations 
to eliminate all artificial interactions between the peri- 
odically repeated supercells due to the dipole moment of 
the slabs. For more details on convergence parameters, 
the construction of appropriate supercell as well as the 
calculated bulk and clean surface structures of ZnO we 
refer to Ref. [jj, where the same computational settings as 
in the present study were used. 

All surfaces were modeled by periodically repeated 
slabs. Very thick slabs consisting of 8 Zn-0 double- 
layers were used to reduce the residual internal electric 
fieldi To represent different surface structures (1x2), 
(1x3) and (2x2) surface unit cells with different combi- 
nations of O vacancies and H adatoms were considered. 
All atomic configurations were fully relaxed by minimiz- 
ing the atomic forces. 

In Table [fl we compare the computed binding ener- 
gies of different bulk and molecular reference structures 
with experimental heat of formations. While the calcu- 
lated binding energies for isolated H2 molecules and bulk 
Zn agree quite well with experiment, there is a notice- 
able error of 1.2 eV in the binding energy of the free O2 
molecule. This is a well known deficiency of DFTi22iSS 
The overestimation of the O2 binding energy is also re- 
flected in a formation energy for ZnO which is to low 
by 0.6 eV. Such deviations would seriously influence our 
analysis of the surface energies, Eq. J2J), and would al- 
ter the allowed range for the oxygen chemical potential, 
Eq. Therefore, to circumvent errors introduced by 
a poor description of the O2 molecule, we have applied 
the following procedure: we take the experimental value 
for the formation energy Ef of ZnO from Table [Q and 
we use Eq. (JSJ to replace the total energy of O2 by Ef 
and the total energies of Zn and ZnO. Ez n and EznQ are 
both bulk quantities which are typically more accurately 
described in DFT than molecular energies. 
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FIG. 1: Contour of the Fermi energy level for the partially 
occupied (0001 )-0 surface bands plotted within the surface 
Brillouin zone (shown to scale). Two 0-2p-bands cross the 
Fermi level (thick solid lines). The unoccupied regions of the 
two bands are indicated by light gray and dark gray areas, 
respectively, bi and b2 are the reciprocal lattice vectors, and 
the surface Brillouin zone is shown by dashed lines. 



III. RESULTS AND DISCUSSION 
A. Surface Distortions 

First we explore the possibility whether the ideal, 
truncated-bulk-like (0001)-O surface may lower its en- 
ergy by breaking the symmetry of the surface layers, 
thereby adopting a distorted surface structure according 
to mechanism (lb). A tendency for symmetry breaking 
reconstructions is often indicated by a strong nesting of 
the Fermi surface. In Fig. ^ we have plotted the two- 
dimensional Fermi surface formed by the partially occu- 
pied 0-2p bands. The plot represents a cut through the 
Brillouin zone of our supercell including only k-vectors 
with a zero component perpendicular to the surface. Fig- 
ure n reveals that actually two surface bands cross the 
Fermi level. This is well known and in full agreement 
with band structure plots presented in Refs. OJ and 0. 
The corresponding wave functions are strongly localized 
at the oxygen atoms of the terminating surface layers and 
are mainly formed by the two 0-2p orbitals parallel to 
the surface. By integrating the occupied and unoccupied 
areas of the Brillouin-zone we find that indeed roughly 
1/2 electron per surface atom is missing to fill the two 
surface bands. However, both Fermi contours are almost 
spherical and only a weak nesting is present. 

As a second test, we did several calculations in which 
we randomly displaced the surface atoms in the top 
atomic layer of our slabs and started an atomic relax- 
ation. Different slabs with (1x2), (1x3) and (2x2) sur- 
face unit cells were used but in all cases the surfaces 
relaxed back toward a symmetric structure with a 3-fold 
symmetry. The surface energy was always higher than in 
the fully symmetric state, so that no hint for a tendency 
toward symmetry breaking reconstructions was found. 
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FIG. 2: Surface free energy A7 of the polar O- 
terminated (000l)-O surfaces with different concentrations of 
O-vacancies cv as function of the oxygen chemical potential 
A/^o- In the top :r-axis, the chemical potential Afio(T,p) has 
been translated into a pressure scale for the fixed temperature 
of T=800 K. Vertical dashed lines indicate the allowed range 
of A/10: For higher chemical potentials O2 will condensate on 
the surface and for lower values of A/10 metallic bulk Zn can 
form. 



B. Oxygen Vacancies 

In the next step we investigate whether the (0001)-O 
surface is stabilized by creating oxygen vacancies. From 
slabs with (lxl), (1x2), (1x3) and (2x2) surface unit 
cells we have removed one O-atom, thereby creating va- 
cancy concentrations cy of 1, 1/2, 1/3 and 1/4. In Fig. [3 
we have plotted the change of the surface energy A7 
of the four defect structures relative to the defect-free 
surface as a function of the oxygen chemical potential 
A/io- As to be expected from mechanism (II) the de- 
fective surface with 1/4 of the O atoms missing is the 
most stable surface structure over a wide range of chem- 
ical potentials. Translating the chemical potential into 
temperature and pressure conditions (assuming that the 
surface is in equilibrium with an O2 gas phase, see upper 
x-axis in Fig. we see that this will be the most sta- 
ble surface at typical UHV-conditions of surface science 
experiments. However, at oxygen rich conditions (high 
pressure and low temperature), exceeding a chemical po- 
tential of — 1.58eV, the ideal, defect-free surface becomes 
the most stable structure. The other surfaces with 1, 1/2 
and 1/3 vacancy concentrations are higher in energy for 
all chemical potentials and will therefore not be present 
in thermodynamic equilibrium. In particular, it is very 
unlikely that the (1x3) structure observed in the HAS 
experiment 1 * corresponds to a simple missing-row struc- 
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c v 1 1/2 1/3 l/4 (a) l/4 (b) 

Ey [eV] +3.24 +2.88 +2.51 +1.80 +1.58 

TABLE II: Calculated vacancy formation energies Ey per 
O atom for removing oxygen from the ideal surface, forming 
O2 molecules and a surface structure with a O vacancy con- 
centration of cy. (a) for a 6-fold (2x2) arrangement of the 
O vacancies with a separation 2a, (b) for a rectangular O va- 
cancy pattern with distances of 2a and y/3a, a being the ZnO 
lattice constant. 

ture with every third O atom removed from the surface. 

At this point we should emphasize that plots like Fig. [3 
strictly only allow to rule out surface structures which 
are higher in energy than other surface models. Since we 
can only perform calculations for a limited set of surface 
models, it is always possible that a not yet considered 
structure with a lower energy exists. For example, since 
we use periodically repeated surface unit cells of a specific 
size in our DFT approach, all our defect structures are 
perfectly ordered. It is however very well possible, that 
disordered or even incommensurate structures might lead 
to a lower energy. Additionally, there are hints that is- 
land and pit structures like the ones observed for the 
Zn-terminated surface^ may also for the O-terminated 
surface be lower in energy than the ideal surface and the 
surface with 1 /4 vacancy concentration considered in our 
studyi^ 

In Table ITT1 we have listed O vacancy formation ener- 
gies Ey which we have defined in the present context as 
the energy difference between the ideal surface and a va- 
cancy surface structure of concentration cy plus 1/2 Eq 2 , 
i.e. Ey is defined with respect to the total energy Eo 2 
of free oxygen molecules, and not, as usually done, with 
respect to bulk ZnO. Ey depends strongly on the va- 
cancy concentration, indicating a strong interaction be- 
tween the vacancies. Up to cv=l/4 the energy cost to re- 
move O-atoms is modest. This is not surprising since up 
to a vacancy concentration of 1/4 the removal of O-atoms 
supports the charge compensation of the O-terminated 
surface and will result in a better filling of the partially 
occupied 0-2p band. For higher vacancy concentrations, 
however, we start to overcompensate the charge transfer 
which stabilizes the surface. The 0-2p band is full now, 
and we have to start to fill Zn-4s-states in the conduc- 
tion band. Therefore the energy cost to remove more 
O-atoms increases rapidly. 



(a) Ideal (0001)-O surface: 

Site: 'on-top' 'hep-hollow' 'fcc-hollow' 

AE [eV] +3.16 0.0 +0.05 

(b) Hydrogen covered (000l)-O surface: 

Site: 'on-top' 'hep-hollow' 'fcc-hollow' 

AE [cV] +1.78 0.0 -0.02 

TABLE III: Relative stability of the high-symmetry ad- 
sorption sites for (a) the surface O atoms of the ideal ex- 
terminated surface and (b) the OH-groups of the H satu- 
rated surface for a full monolayer coverage. The AE are the 
calculated energy differences per O-atom/OH-group for mov- 
ing the topmost O/OH-surface layer from the regular lattice 
position of the wurtzite structure ('hcp-hollow-site') to the 
'on-top' and 'fcc-hollow' position. 

for different surface models with H adatoms, we consider 
the possibility that the preferred adsorption site of these 
OH-groups is no longer the regular lattice site of the O 
ions. Three different high- symmetry adsorption sites are 
conceivable above the underlying Zn layer: an 'on-top' 
position, a 'hep-hollow site', which is the regular lattice 
position for the O ions in the wurtzite structure, and a 
'fcc-hollow site'. 

First we consider the clean surface without adsorbed 
hydrogen. Then we see from Table IIIII that indeed the 
'hep-hollow' sites are the most stable positions for the 
O surface layer. However, moving the whole layer to 
'fcc-hollow' sites costs only a small amount of energy. 
Turning to the hydroxylated surface with a full mono- 
layer coverage of hydrogen we find that the OH-groups 
now prefer the 'fcc-hollow site'. So by gradually adding 
hydrogen, the regular lattice site of the surface O ions be- 
comes unstable relative to the 'fcc-hollow site'. However, 
the energy difference between 'fee-' and 'hep-hollow sites' 
is so small that we have neglected this effect in all further 
calculations and have only considered 'hep-hollow sites' 
for surface oxygen atoms and OH-groups. 

In the next step we construct different surface models 
of hydrogen covered (0001)-O surfaces using a similar 
procedure as in Sec. ImBl We take slabs with (lxl), 
(1x2), (1x3) and (2x2) surface unit cells, and we add 
different amounts of hydrogen to create hydrogen cover- 
ages of ch of 1/4, 1/3, 1/2, 2/3, 3/4 and 1 monolayer. 



C. Hydrogen Adsorption 

We turn now to a situation in which the (OOOl)-O sur- 
face is in equilibrium with a H2 gas phase. H2 molecules 
can dissociate, and hydrogen atoms may adsorb at the 
surface thereby forming OH-groups with the surface oxy- 
gen ions. Before we start to calculate the surface energy 



c H 1 3/4 2/3 1/2 1/3 1/4 

E h [eV] -0.71 -1.10 -1.25 -1.90 -2.12 -2.20 

TABLE IV: Calculated binding energies per H atom for 
dissociating H2 molecules and forming hydrogen layers of cov- 
erage ch. 



7 



The calculated surface energy changes A7 relative to the 
clean (OOOl)-O surface are plotted in Fig. as a func- 
tion of the hydrogen chemical potential A^h- A very 
similar behavior as in Section IIII Bl arises: At H-rich 
conditions the structure with a 1/2 monolayer hydrogen 
coverage is the most stable surface, in agreement with 
mechanism (III). On the other hand, at H-poor condi- 
tions the clean surface without hydrogen becomes the 
most stable structure. As a new feature we find that be- 
tween the two limiting cases the surfaces with 1/3 and 
1/4 monolayer coverage are slightly more stable for small 
intervals of the chemical potentials,— Thus, by lowering 
the chemical potential A/in from H-rich to H-poor con- 
ditions it is possible to gradually reduce the hydrogen 
coverage from 1/2 monolayer to 1/3, 1/4 and zero cover- 
age. Translating the chemical potential into temperature 
and pressure conditions we find that we will start to re- 
move H at UHV-pressures at a temperature of roughly 
750 K, which is in reasonable agreement with the exper- 
imental observation^ The other surface structures with 
hydrogen coverages larger than 1/2 are always higher in 
energy and will be unstable for all temperatures and hy- 
drogen partial pressures. In particular a surface with a 
full monolayer of hydrogen as postulated from the results 
of the HAS experiments^ is not likely to exist in thermo- 
dynamic equilibrium. However, from the intensity of the 
He-specular peak it was deduced that only about 0.1 % 
of the (0001)-O surface consists of flat terraces with di- 
ameters exceeding 50 A which contribute to the (lxl) 
HAS signal i& Therefore, the H covered surface with a 
(lxl) HAS diffraction pattern may well be a minority 
phase which is formed under suitable kinetic conditions. 

In Table llVl we have summarized the H binding ener- 
gies .Eb per atom which we have calculated as energy dif- 
ference between the ideal (0001)-O surface plus 1/2 E# 2 
and surface structures with H coverages of ch- We see 
that it becomes rapidly unfavorable to adsorb more hy- 
drogen as soon as the concentration for ideal charge 
compensation of the polar surface of 1/2 monolayer is 
reached. The reason for this behavior is the same as in 
the case of the oxygen vacancies: Up to 1/2 monolayer of 
hydrogen we fill the partially occupied 0-2p-band, be- 
yond 1/2 monolayer the 0-2p-band is completely filled 
and we have to populate the conduction band. The de- 
crease in E\, from 1/2 to 1/4 monolayer coverage indi- 
cates that also at low coverages a weak repulsive inter- 
action between the hydrogen atoms remains. This is the 
reason why also the low-coverage structures appear in the 
surface phase diagram. If we extrapolate toward zero 
coverage of the surface we can estimate an initial bind- 
ing energy of about 2.3 eV per H atom for the dissociative 
adsorption of H 2 . 



D. Water Dissociation 

As evident from Fig. [31 hydrogen adsorption plays a 
major role for the stabilization of the polar (0001)-O 
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FIG. 3: Surface free energy A7 of the polar O-terminated 
(000l)-O surfaces with different coverages of hydrogen ch as 
function of the hydrogen chemical potential A^ih- In the top 
:r-axis, the chemical potential A[m{T, p) has been translated 
into a pressure scale for the fixed temperature of T=800 K. 
The vertical dash line indicates the upper bound for A/_tH. 

surface and for almost all conceivable experimental con- 
ditions hydrogen will be present on the surface. But until 
now we have only considered molecular hydrogen as the 
reservoir for the surface hydrogen. However, in many 
chemical reactions and catalytic processes also water is 
present. Therefore we will briefly explore if also water 
can act as a source for surface hydrogen in the presence 
of the (0001)-Zn face. 

In a DFT study using the hybrid B3LYP functional 
Wander and Harrison-? found that dissociating water and 
forming a full H and OH monolayer on the O- and Zn - 
terminated surface, respectively, is energetically unfavor- 
able by 0.1 eV. As adsorption sites for the OH-groups 
they assumed the Zn 'on-top' positions which would be 
the next lattice sites for O ions if the crystal is ex- 
tended. However, on the polar surfaces two more high- 
symmetry adsorption sites exist: the 'hep-hollow site' 
position above atoms in the second surface layer and a 
'fcc-hollow' site with no atoms beneath. 

Using the same adsorption geometry as Wander and 
Harrison we also find that dissociating water is energet- 
ically unfavorable with a slightly larger energy cost of 
0.3 eV. However, as shown in Table Ivl the configuration 
with the OH groups adsorbed at the 'fcc-hollow sites' in- 
stead of the 'on-top' positions is much lower in energy. 
Considering the correct adsorption positions for the OH 
groups we now find that even for the thermodynamically 
unstable monolayer coverages the dissociation of water is 
energetically preferred by about 0.4 eV. Taking only 1/2 
monolayer coverages into account, this energy gain will 
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Site: 'on-top' 'hep-hollow' 'fcc-hollow' 

AE [eV] 0.0 -0.04 -0.72 

TABLE V: Relative stability of the different OH adsorp- 
tion sites on the polar Zn-terminated surface, calculated for 
a monolayer coverage. 



be significantly larger. 



E. Surface Phase Diagram 

Finally we combine the results of the previous subsec- 
tions and assume that the polar (000l)-O surface is now 
simultaneously in equilibrium with an O2 and a H 2 gas 
phase. In addition to the surface models described in 
Sec. 1111 Bl and Sec. IIII CI we have furthermore considered 
various mixed structures of O vacancies and adsorbed H 
atoms in the (1x2), (1x3) and (2x2) surface unit cells. 

The surface free energy now depends on two chemi- 
cal potentials A/iq and A/in- The graphs of Fig. [21 and 
13 therefore have to be extended to a 3-dimensional plot. 
Such a diagram would be rather complex and hard to fol- 
low. The most important information contained in the 
plot of the surface free energies is which of the surface 
models has the lowest surface energy for a given combina- 
tion of chemical potentials A/io and A/zh- This informa- 
tion is better visualized if we project the 3-dimensional 
diagram onto the (A/io,A/ih) plane and only mark the 
regions for which a certain surface structure is the most 
stable one. The result is a phase diagram of the (0001)-O 
surface which is shown in Fig. 01 

The surface phase diagram in Fig. 0] summarizes in a 
condensed fashion the essential results of our study. This 
phase diagram is dominated by two structures: a surface 
with 1/2 monolayer of adsorbed H and a hydrogen- free 
surface with 1/4 of the oxygen atoms removed. These are 
the two scenarios denoted as mechanism (II) and (III) in 
the Introduction, indicating that filling the 0-2p-bands 
is a very important mechanism to stabilize the (0001)-O 
surface. Next to these two phases we find two structures 
in which H is gradually removed from the surface and 
only at very H-poor conditions and when plenty of oxy- 
gen is available, the ideal O-terminated surface stabilized 
by mechanism (I) becomes the most stable structure. 

None of the additionally considered mixed structures 
which simultaneously contain O vacancies and H adatoms 
appears in the phase diagram. This is not very surprising 
since for the given sizes of the surface unit cells no com- 
bination of O vacancies and H adatoms exists that leads 
to fully occupied 0-2p-bands. However, it is very well 
conceivable that for larger surface unit cells mixed struc- 
tures with, for example, an O vacancy concentration of 
1 /8 and a H coverage of 1/4, become more stable which 
would then appear as new phases between the H covered 
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FIG. 4: Phase diagram of the polar O-terminated (0001) sur- 
face in equilibrium with H and O particle reservoirs control- 
ling the chemical potentials Afiu and A/io, based on selected 
superstructures as explained in the text. The lowest-energy 
surface structures are labeled by the concentrations of oxy- 
gen vacancies cy and hydrogen adatoms ch . The upper right 
area indicates conditions under which H2O condensates on 
the surface. 



and the O vacancy structures. 

Relating the chemical potentials to temperature con- 
ditions and partial pressures of the gas phase shows, see 
Fig. 01 that for almost all realistic experimental condi- 
tions hydrogen will be present at the (0001)-O surface. 
Even at UHV-conditions with a low hydrogen partial 
pressure one has to go to rather high temperatures to 
fully remove all hydrogen. In this case, a surface struc- 
ture with O vacancies will become the most stable one. 
In order to stabilize the ideal O-terminated surface an 
oxygen atmosphere with an extremely low content of hy- 
drogen (and also water vapor) is necessary, which is ba- 
sically not achievable in experiment. 

In Table IVII we have summarized the surface relax- 
ations for the three most important surface structures 
appearing in the surface phase diagram, Fig. 01 For the 
extended surface structures with H adatoms and O va- 
cancies we have averaged in each atomic plane parallel to 
the surface the atomic displacements, and we define the 
interlayer distances d as the separation of the averaged 
atomic positions. Depending on the surface structure 
and the charge compensation process, very different sur- 
face relaxations occur. The largest relaxations are found 
for the clean, defect-free surface termination with a com- 
pression of the first double-layer distance of almost 50 % 
and also a significant contraction of the second and sub- 
sequent double-layer spacings. This is in agreement with 
other previous ab-initio studies iL&iii After filling the 
partially occupied surface bands by adsorbing 1/2 mono- 
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d ideal surface H covered O vacancies 



this surface phase diagram we predict that hydrogen is 
adsorbed at the (OOOl)-O surface for a wide range of 
temperatures and H 2 partial pressures, including UHV- 
conditions. This is in agreement with the recent observa- 
tions of a HAS experiment^ and was also confirmed in a 
study of the CO adsorption on the polar ZnO surfaces. 19 

We find a H binding energy of roughly 2.3 eV per 
atom if molecular hydrogen dissociates and adsorbs at 
the clean O-terminated surface. Furthermore we predict 
that in situations where both polar surface terminations 
are present (for example for powder samples) also the 
dissociative adsorption of water with H and OH-groups 
being adsorbed at the O- and Zn-terminated surface, re- 
spectively, is energetically preferable. However, as soon 
as a coverage of 1/2 monolayer of hydrogen is reached, the 
energy gain of adsorbing more hydrogen on the (0001)- 
O surface drops very rapidly with increasing hydrogen 
coverage. Therefore no stable phases with more than 
1/2 monolayer H coverage appear in the surface phase 
diagram, Fig. 0] In particular, a structure with a full 
monolayer of H as predicted in Ref. [18| is not very likely 
to exist globally in thermodynamic equilibrium (which 
does not exclude a kinetic or local stabilization). 

Going to low hydrogen partial pressures and higher 
temperatures it is possible to gradually remove the hy- 
drogen from the surface and to form stable phases with 
less than 1/2 monolayer coverage of hydrogen. If all hy- 
drogen is removed, ox yge n vacancies will be created as 
was speculated in Ref. [3 However, we find that a sur- 
face with a vacancy concentration of 1/4 is much more 
stable than a missing-row structure where 1 /3 of the oxy- 
gens has been removed. Therefore, based on the limited 
set of surface structures taken into consideration in our 
study, we currently do not understand the (1x3) struc- 
ture observed in Ref. Ill 

Finally, at higher oxygen partial pressures the O vacan- 
cies will be filled and the clean, defect-free O-terminated 
surface becomes the most stable structure. However, the 
hydrogen partial pressure has to be very low so that we 
consider it very unlikely that a clean, defect-free (0001)- 
O surface can be observed in experiment. 



V. ACKNOWLEDGMENTS 

The author would like to thank Dominik Marx, 
Christof Woll, Georg Kresse, and Ulrike Diebold for fruit- 
ful discussions. The work was supported by SFB 558 and 
FCI. 



H-Oi 0.1825 

Oi-Zn 2 0.0628 (-48%) 0.1207 (0.0%) 0.1151 (-4.6%) 

Zn 2 -0 3 0.3985 (+5.1%) 0.3779 (-0.4 %) 0.3767 (-0.7%) 

3 -Zn 4 0.1077 (-11%) 0.1246 (+3.2%) 0.1238 (+2.6%) 

Zn 4 -0 5 0.3813 (+0.5%) 0.3773 (-0.5 %) 0.3772 (-0.6%) 



TABLE VI: Summary of the averaged interlayer distances d 
(given in fractions of the theoretical bulk lattice constant c) 
and the relative changes with respect to the bulk values (in 
parentheses) for three different surface structures: the ideal, 
defect-free (0001)-O surface, the H covered surface with 1/2 
monolayer of hydrogen and the surface with a vacancy concen- 
tration of 1/4. The subscripts refer to surface layers numbered 
from the surface plane. The theoretical PBE bulk values for 
the interlayer distances are do-z n =0.1208 c (in bilayer) and 
dzn-o=0.3792 c (between bilayers) and the lattice constant 
was calculated to be c=5.291 A (see Ref. 0). 



layer of H or by removing 1/4 of the O ions, however, 
the surface layers relax back to almost truncated-bulk- 
like positions. Thus, for surfaces with a lower H adatom 
or O vacancy concentration, surface relaxations between 
the two extremes are conceivable. This may explain why 
experimentally very different results for the surface re- 
laxations were found. GIXD measurement a 8 ! 11 have pre- 
dicted an inward relaxation of the upper O-layer with 
a contraction of the first double-layer distance of 40 % 
and 20 %, respectively, whereas from LEEDiS and LEIS 12 
experiments it was concluded that the first double-layer 
spacing is close to its bulk value. 



IV. SUMMARY AND CONCLUSIONS 

By combining first-principles density-functional cal- 
culations with a thermodynamic formalism we have 
determined lowest-energy structures of the polar O- 
terminated (0001)-O surface of ZnO in thermal equilib- 
rium with an O2 and H2 gas phase. This scheme allows us 
to extend our zero-temperature and zero-pressure DFT 
results to more realistic temperature and pressure condi- 
tions which are usually applied in surface science exper- 
iments or in catalytic processes like the methanol syn- 
thesis, and thus to bridge computationally the 'pressure 
gap'. 

The essential result of this approach is a phase diagram 
of the (0001)-O surface which is plotted in Fig. From 
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